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PREFACE 



Two feedback control schemes which maximize a terminal quantity while 
satisfying specified terminal conditions are discussed and compared. The 
schemes are based on a linear perturbation from a nominal non-optimal path 
which does not, in general, satisfy the terminal conditions. The methods 
have been programmed in the ALGOL computer language for evaluation and the 
programs are included in the appendices. 
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I . INTRODUCTION 



In recent years several authors have treated the problem of deter- 
mining optimum control programs for nonlinear systems with terminal 
constraints. These problems arise in the design of control systems and 
development of guidance laws where it is desired to determine, out of all 
possible time histories of the control variables, the one control history 
that maximizes (or minimizes) one terminal quantity or cost function while 
simultaneously yielding specified values of certain other terminal 
quantities . 

The steepest-ascent or gradient method developed by Kelley\ Bryson 
2 

and Denham , which is a systematic and rapid numerical procedure, has 
proved to be successful in solving this class of problems. Improvement 

in the convergence time of the iterative process involved has been achieved 

3 k 

by Rosenbaum by a method based on the earlier work of Bryson and Denham . 

Another successful method developed by Breakwell, et.al.^, and modified 

g 

by Bullock is a second variation method in the Calculus of Variations. 

The principal objectives of this thesis are to develop a simpler 
steepest-ascent program which will be understandable to the control 
engineer without a background in the Calculus of Variations and to compare 
the results and speed of convergence with the method developed by Bullock. 

In this way it is expected that the steepest-ascent program will prove to 
be a useful instrument in education and research, while at the same time 
through the comparison, illustrate the advantages and disadvantages of the 
two approaches to the problem. 

The method used in developing the steepest-ascent program is 

4 

essentially a variation of the methods of Bryson and Denham , and 
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Rosenbaum ? hence it is restricted to problems in which only the de\iations 
in the control variables and adjustable parameters are considered in the 
performance index. It is also restricted to problems in which the pay-off 
quantity is a function of the terminal value of the states. This varia- 
tion includes a terminal- error control scheme which maintains a bound on 
the terminal constraint errors , hence reducing the total number of itera- 
tions required to converge to the optimum since larger deviations from the 
nominal trajectory can be tolerated while still meeting the desired terminal 
conditions . 

A numerical example is given of a rocket ascent trajectory into a 
circular orbit of maximum altitude. Provision is made for a two-stage 
rocket with optimization of the inter-stage coast duration. 
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II. STATEMENT OF THE PROBLEM 



Given a system which can he described by a set of non-linear (or 
linear) first order, ordinary differential equations, determine a control 
history u(t), in the interval S t S T, to maximize 

P = x^T), (1) 



subject to constraints 

t = *[x(T)] = 0, (2) 

= f[x(t),u(t),t], (3) 



t Q , T, and x(t Q ) given; 



where 



u(t) 



\(t) 



u (t) 
m v 7 



an m x 1 matrix of 
control variables > 



oo 



(5) 



x(t) 



x i(t) 
X (t) 

n 



an n x 1 matrix of state 
variables , 



( 6 ) 



4 ' 






a q x 1 matrix of terminal 
constraint functions , each 
of which is a known function 
of x(t) , 



(7) 



f 



r 

n 



an n x 1 matrix of known 
functions of x(t)> u(t), 
and t, i.e., the system 
equations; 



( 8 ) 
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( 9 ) 



f> is the pay-off quantity and is one of the 
states, namely x^(t); 

T is the terminal value of the independent variable. 

On each iteration it is desired to minimize the mean value of a 
positive definite quadratic form in the control variable deviations: 



where the superscript ( ') indicates the transpose of a matrix, and 6 u(t) 
are small deviations in the control history from a nominal non-optimum 
trajectory. 




T 



( 10 ) 



h 



III. STEEPEST -ASCENT METHOD OF SOLUTION 



A. BACKGROUND 

i 

Bryson and Denham, in Ref. considered terminal control of non-linear 
(and linear) systems for minimum mean values of a positive definite quad- 
ratic form in the control variable deviations. That is, it was assumed 
that a nominal control history had been determined which caused the vehicle 
to arrive at the terminal point with desired values of certain specified 
terminal conditions. Small deviations from this nominal trajectory were 
considered which might be caused by disturbances, inaccuracies in the ■ 
data, inaccuracies in the control system, etc. The problem was to determine 
small deviations from the nominal control so that the terminal constraints 
would be satisfied in spite of the disturbances. 

In the present paper the nominal trajectory is determined by guessing 
a reasonable control variable program. For example, in a rocket trajectory 
problem one might choose an initial launch angle and a gravity turn with 
zero thrust angle throughout as is done in the numerical example. Further- 
more, it is desired not only to determine control deviations which result 
in meeting the terminal constraints, but also to maximize the terminal 
value of one of the states while minimizing a performance index. This 

optimization scheme is a variation of the so-called Lambda Matrix Control 

A 3 

feedback method described in Ref. / and the convergence method of Ref. j?. 

B. DERIVATION OF EQUATIONS 

The optimum programming problem can be solved systematically and 
rapidly on a high speed digital computer using the steepest-ascent 
technique. As stated in (A), this technique -starts by guessing a nominal 
control variable program, u*(t), and solving the set of differential 
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equations (3) with initial conditions (4), to determine a nominal tra- 
jectory. This trajectory, in general, will neither maximize f) nor will it 
satisfy the terminal constraints (2). 

Consider small perturbations in the control variables, &u(t), about 
the nominal, where 

5u(t ) = u(t) - u*(t) (n) 

(The superscript (*) indicates terms evaluated along the nominal trajectory.) 
These perturbations will cause perturbations in the state variables, &x(t), 
where 

5x(t) = x(t) - x*(t) (12) 



Substituting these relations into the system differential equations (3) and 
expanding f in a Taylor series about the nominal the result is, to first 
order 



ft (6x) = F(t)6x(t) + D(t)6u(t) 



(13) 



where 



and 




an n x n matrix of 
partial derivatives 
evaluated along the 
nominal trajectory; 



D(t) 



-(0 



an n x m matrix of 
partial derivatives 
evaluated along the 
nominal trajectory 



To determine the effects upon the terminal conditions 0 and \|/ we 



introduce the linear differential equations adjoint to (12) defined as 



— (T,t) = -®(T,t)F(t) 



6 



(14) 



where 0 is an n x n fundamental or state transistion matrix whose elements 



give the sensitivities of the terminal states to perturbations, Sx(t), 
along the trajectory. (See Ref. 7)* Initial conditions for these equations 
are specified at the terminal time, i.e., 

$(T,T) = I, the identity matrix (15) 

hence numerical integration proceeds backward in time. 

The solution to (l4) provides a solution to the linear perturbation 
equations (12) at the terminal point: 



6x(T) = o(T,t)dx(t) 




®(T, T)D(r)Su(T)dT 



( 16 ) 



Since f> and i); depend on the terminal values of the states, small 
deviations, aft and d\|i may be calculated from the solution to (l6). Con- 
versely, if 5 x(t) is known by specifying values of d^ and di|r, the corres- 
ponding control history may be calculated, which is what is done. Re- 
writing ( 16 ) 




(17) 



In order to minimize the performance index subject to the constraints 
(16), the method of Lagrange multipliers is employed. Multiplying (l6) by 
a matrix of Lagrange multipliers 




(18) 
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(where (l6) is written to include only those states which appear in f> and 
\|» 3 hence q+1 equations) and adjoining the result to (10) 



C - v 'Sx(T) - v'®(T,t)Bx(t) +J 8u’(t)5u(t) - v , <D(T 5 t)d(t)Bu(t)J dr 



( 19 ) 



Completing the square 



C = v'[&x(t) - $(T,t)Sx(t)] +f ft (6 u(t)D ' (t)3> ' (Tjt)v) ' (6u(t) -D ' (t)$ '(T,t)v) ■ 

A L 

- ^ v '$( t 3 t)d(t)D ' (t)$ ' (T 3 r)vjdr (.20) 



To minimize C subject to the control variations choose 

5u(t) = D '(t)$ '(Tjt)v (21) 

Substituting this relation in (17) in order to find v 



r T 

Bx(t) - C'(T,t)Sx(t) -J $(T,T)D(T)D'(T)$-(T 3 T)vdT = 



( 22 ) 



Define the controllability matrix 



= + L 



$(T 3 t)d(t)D'(t)$ '(T 3 t) d T 



(23) 



Note: This integration may be performed simultaneously with (l4) by 

numerically solving 



|j: = -0(T 3 t)D(t)D'(t)5>'(T 3 t) 



(24) 



with initial conditions 

J(T) = 0 (25) 

Equation (22) may be written 
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5x(T) - C>(T,t)5x(t) = Jv 



(26) 



Solving for v 

v = J -1 [Sx(t) - 4>(T,t)6x(t)] (27) 

Substituting this relation in (2l) 

Bu(t) = D , (t)0' , (T,t)j' 1 [5x(T) - C>(T,t)Sx(t)] (28) 

which is the perturbation in the control history which satisfies the con- 
straint (17) while minimizing the performance index (lO). 

As stated earlier, 5 x(t) is determined from the specified values of 
dcp and d\[r where 



dcp = 6x. (T ) 


(29) 


d* = *[8x(T)] 


(30) 



As yet, nothing has been said as to how one chooses the desired pay- 
off improvement, dcp, or the desired improvement in meeting the terminal 
constraints, d\|/. The latter is normally chosen as 

d\[r = -i|f (31) 

that is, the negative of the total error on any iteration is chosen as 
the desired correction specified on the following iteration. The problem 
of specifying dcp is more complex and is the subject of the next section. 

It is worthy to note at this point that the Lagrange multipliers, 
which are error feedback terms b need not be computed at every point on 
the trajectory. Sufficient accuracy can be obtained in computing the 
control deviations (21) by solving (27) at discrete intervals and using 
the result until the next "sampling time". This reduces the number of 

9 



times the controllability matrix must be inverted per iteration and materi- 
ally improves the running time of the program. Experimentation will re- 
veal how large the sampling interval can be made. Since the controllability 
matrix is singular at the terminal time, T, new values of the Lagrange 
multipliers should not be computed too close to the end of the trajectory. 



C. METHOD OF SPECIFYING THE IMPROVEMENT IN PAY-OFF 

In general, one does not know how far from the optimum a given nominal 
trajectory will be. It is, therefore, difficult to guess how much pay-off 
improvement to specify initially. However, it is possible to compute a 
value of dcp which will result in a trajectory that satisfies the terminal 
constraints. This is done as follows: 

The changes in the control variables required to meet the terminal 
constraints with the pay-off unconstrained are given by 



where 




1-1 



A A f dT 



d\|/ 



A = $(T,t)D(t), without row 1 or column 1. 



(32) 



This is the same as the basic control equation (28) with ^(t^) equal to 
zero and 5 x(t) containing only the terminal constraint terms. 

The change in pay-off, dcp, that will be produced by a given change 
in the control variables is, from (l6) 



dcp = 




( 33 ) 



where is the first row of the A matrix. Substituting (33) in (32) 
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where is the first row of the A matrix. Substituting (33) in (32) 



dcp = 




( 34 ) 



Equation (34) gives the change in pay-off associated with adjusting 
the control in order to meet the terminal constraints. This value is used 
on the first iteration. 

Equation (34) is also used to compute a value of the pay-off cor- 
rected for terminal error s 5 i.e. 5 cp + dcp is the value the pay-off would 
have achieved had the terminal errors been zero. Hence 5 one can determine 
whether an improvement in the pay-off was actually achieved or if an 
apparent improvement was a result of larger terminal constraint errors. 

On subsequent iterations, one of three methods is used to compute dcp. 
A value equal to 25 per cent of the nominal value of cp is computed and 
stored. This quantity is called dcp*-* and is used in method (2) helow. It 
is a fairly arbitrary choice hut should he made as large as seems reason- 
able. The program will automatically adjust it if it is too large. 

Method 1. Choose dcp to satisfy the terminal constraints 

with the pay-off unconstrained as described above. 

Method 2. If |d^| < e, where s is chosen as reasonable tolerance on 
the terminal constraints then 

acp=^ 

2 1 

where i is a count of the number of times method 
(3) has failed. This has the effect of halving 
the improvement specified each time a run is 
unsuccessful in improving the terminal errors or 
the corrected value of the pay-off. The program 
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terminates when, while executing this method, 
dcp becomes less than a pre-set number. A 
final run is then made using method (l). 

Method 3. If |d\|f| >e the following questions are asked: 
Were the errors on the current iteration 
smaller than on the preceeding iteration? 

Was there an improvement in the corrected value 
of the pay-off? If the answer to either 
question is no, the run is considered •unsuc- 
cessful, the control history is replaced by the 
previous control history, and method (l) is 
used. If the answer to either . question is yes, 
dcp is set equal to zero and an attempt is made 
to satisfy |d\p| < e. If this test fails a 
second time, method (l) is used. 



D. TWO-STAGE ROCKET TRAJECTORY WITR COAST PERIOD OPTIMIZATION 

In many orbit injection applications, such as the Gemini-Titan II 
system, the launch vehicle is made up of two powered stages. It is there- 
fore of interest to consider the effect of an interstage coast phase on 
the maximum altitude obtainable. In this section a method of calculating 
the optimum coast duration is derived. 

The basic equations, (l) through (12), are the same. The linear 
perturbation equations (13) may be written 

[Sx(t)] - F(t)5(t) + D(t)8u(t) + B(t)5c 

where 

an n x 1 matrix of partial 
/ \ _ derivatives of f with 

^ ' \ 5c / 5 respect to the coast 

duration, c. 



( 36 ) 



( 37 ) 
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The solution to (37) is 




T $(T,T)B(T)8cd T (38) 



The performance index becomes 




(39) 



^ •'t ^ 

where A is simply a weighting factor. 

Before attempting to minimize (39) subject to the constraint (38), 



change in the terminal values of the states due to a change in the coast 
duration. Since c is an adjustable parameter which dees not appear ex- 
plicitly in f, the partial derivatives, cannot be evaluated directly. 
However the desired term may be calculated by the following method: 
Define 

t^ = Stage I burnout time 
t^ = Stage II ignition time 

hence 



a method must be derived to evaluate the last term in (38) which is the 




(40) 



Since 5c is a small time increment we may write 




-x(t2)+Ax 

x(t 2 ) 



= (x(t 2 ) + Ax) - x(tg); 
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consequently we may make the approximation 



x(tg + 5c) - x(t 2 ) = (—] 5c (4l) 

Now define x^ as the states evaluated with the thrust off (uncontrolled), 
and x c as the states evaluated with the thrust on (controlled). From (4l) 

x u ( t 2 + Sc) = x(t 2 ) + x u (t 2 )Sc 

x c (tg + 5c) = x(tg) + x c (t 2 )6c 

where (*) indicates differentiation with respect to time. Subtracting the 
above expressions we have 

x (t 0 + 6c) - x (t 0 + 6c) = [x (t ) - x (t 0 )]6c 
u' 2 ' c 2 7 u' 2 7 c v 2 7 

Define 

6x c = f x u^2^ " ^ c ^ t 2^ 5c ' (^ 2 ) 

The quantity 6x^ is a perturbation in the states occurring at time 
t = t 2 + 6c 

due to cutting off the thrust for a period 6c. The question remains: 

How does this perturbation propagate to the end of the trajectory? The 
answer is clearly 

6x(t) = <t(T,t 2 + 6c)6x c (43) 

Finally, (38) may be written 

r T 

6x(T) = «>(T,t)6x(t) + / <p(T,T)D(-r)Su(T)dT + $(T,t + 5c)&x (44) 

J t 2 c 

where the last term is zero prior to time t^ *+ Sc. 



Introducing the Lagrange multipliers, v, and adjoining ( bk ) to (39) 



/* T 

C = | J Su ' (t)5u(t)cIt + | ASc 2 + v’6x(T) - 6 ’<f(T,t)&x(t) 



- V 



r T 

J $(T,t)d(t)&u(t) dr - v '0>(T,t 2 + 6c)&x c 



(45) 



As before, completing the square yields the optimum control change 

5u(t) = D t (t)0 ! (t,t)v (46) 

By differentiation the optimum coast change is 

6c = y v ' $(T,t 2 + 6c)[x u (t 2 ) - x c (t 2 )] (47) 

Substituting (46) and (47) into (44) 



5x(l) = 6(T,t)6x(t) + f ^(T,T)D(r)D'(,T)$'(T,T)dT 

+ ^ v '$(T,t 2 + 6c)[x u (t 2 ) - x c (t 2 )][x u (t 2 ) - x c (t 2 )]' 0'(T,t o +6c)i 



Define 



A = &x(t) - $(T,t)6x(t) 
r T 

J =/ $(T,T)D(T)D'(T)©(T > T)dr 

A = fx u (t 2 ) - x c (t 2 )] 

G = $(T,t 2 + 6c )A A’ 0'(T,t 2 + 5c) 



(49) 

(50) 

(51) 

(52) 



Rewriting (48) 



A = 



Jv + j- Gv 



V 




(53) 

(54) 
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Finally 5 substituting in ( 46 ) and (47 ) 5 the control and coast variations are 




( 56 ) 



( 55 ) 



Since the last term in (44) is zero prior to time t^ + 5 , the term 
G/A in (55) and ( 56 ) is also zero prior to that time. 



Equation (5 6 ) should he evaluated at t = tg + 5c, but since 5c is the 
unknown, it is evaluated at t^. This does not introduce an appreciable 



error if 5c is small. 

In some numerical integration procedures the terminal value of the 
independent variable cannot be changed once the integration has begun, 
hence the change in coast time cannot be added immediately. This problem 
is solved by evaluating ( 5 6 ) on each forward integration (except the nominal) 
and, if 5c f 0, re -integrating the latter portion of the trajectory from 
t^ to T. 

E . COMPUTATIONAL PROCEDURES 

As stated earlier, this steepest-ascent technique requires the use of 
a high speed digital computer. The sequence of operations is summarized 
here . 

1. Compute the nominal path by integrating the system differential 
equations (3) with a nominal control history and appropriate initial con- 
ditions and store the time history of the state variables at reasonably 
small intervals. Print out the values of cp and i|r. 

2. Integrate the adjoint differential equations ( 14 ) backward, 
evaluating the partial derivatives on the nominal path by reference to 

the states stored in step (l). Simultaneously integrate the controllability 
matrix equations (24). Store the results at the same interval as the states. 



3. Select desired terminal condition changes , dcp and d^, as ex- 
plained in Sections B and C. 

b. Compute and use the new control history while integrating the 
system differential equations forward. Again print out the values of cp 
and \|/ 5 unless the next step applies. 

5- If the two-stage rocket problem is being solved, compute the new 
coast period in step (4). Transfer the storage locations of the second 
stage control history to correspond with the new coast time. Integrate the 
system equations from t^ to the new terminal time. Print out the values of 
cp and >|r. 

6. Repeat procedures (2) through ( 5 ) until the pay-off improvement 
in step (3) is less than a preset value. At this point, use method (l) 
described in Section C to select dcp and complete step (4) and (5). This 
has the effect of eliminating any remaining errors in the terminal 
constraints . 

7 . Punch cards or store the control history on tape and terminate 
the program. 

Before concluding this section, a few general factors of great 
importance in this type of numerical calculation should be discussed. 

The programmer must exercise great care when working with values of 
type real (or floating point). Often a calculation is made where the 
result is expected to be an integral value such as 4/2 = 2.000 ..., 
however, due to the binary, octal, and decimal conversions which take 
place within the computer, the result may come out 1.9999- • *99* This 
problem occurs when trying to generate array storage indices based on a 
value of the independent variable which is a floating point quantity. 
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IV. SECOND VARIATION METHOD OF SOLUTION 



A. OUTLINE OF THE METHOD 

Bullock^ has derived a feedback control scheme based on the second 
order variational theory in the Calculus of Variations. The method is 
outlined here in sufficient detail to solve the problem stated in II. 
The differential equations to be satisfied are 



X = f(x, u, t) 



(57) 




(58) 




(59) 



and 



b = M’d - N'c 



(60) 



in order to maximize 



9 = <P[x(T), T] 



(61) 



and satisfy the terminal constraints 



♦ = t[x(T), T] 



(62) 



Define the variational Hamiltonian 



H = A'f 



(63) 



The elements of Eqs. (59) and (60) are 



F = f - f H' 1 H 



(6k) 



X u uu ux 




( 65 ) 
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( 66 ) 



S » H - H H ' 1 H' 
XX XU uu u 



c = - f h " 1 H' 
u UU u 



d = H H ' 1 H' 
XU UU u- 



(67) 

(68) 



where the subscripts indicate partial differentiation in the usual sense. 

The initial conditions for (57) are given and the terminal conditions 
for (58), ( 59)5 and ( 60 ) are 

= (q> - V'lr ) (69) 

x x t=T 

where the components of the column vector v are sensitivities of the pay- 
off <p to changes in the terminal constraints \|r; 



m(t) = i - ru ry 1 \|r ( 70 ) 

JL -A. X. 

where I is the identity matrix; 

N(T) = (71) 

b(T) = -i|r x dx)r (72) 

where d\|r = -\|r. 



The perturbations in the control history are given by 



5u(t) = u* - u = -H _1 (H + H Sx + f ’ SA) ( 73 ) 

v J uu v u ux u J ' 

where 

6A = (M 1 ) _1 (N '6x + b). (74) 

In order to minimize a performance index 
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( 75 ) 



C = ~ f 6u'(T)6u(T)dT 

H in the above equations is modified, by adding an arbitrarily large 
negative constant, K, which is reduced in magnitude as the program con- 
verges. This has the effect of constraining the magnitude of 5u. 

Since the terminal conditions on the adjoint equations, A(r) , depend 
initially on the choice of v, a method is given which will improve the 
accuracy of these terminal conditions on subsequent iterations. 

Equation (7^-) is an expression for 6A at any time, t, but since M is 
singular at T, it cannot be evaluated directly. However, if a point (t^) 
is chosen sufficiently far from this singularity, the following equation 
can be integrated from t-^ to T: 

SA = -SSx - F* SA (7 6 ) 

with initial conditions 

SA(t, ) = (M')' 1 (N’ 8 x + b) (77) 

1 t=tl 

The solution to (76) is then added to the current values of ?\ ( T) prior to 
the next backward integration of ( 58 )* 

B. COMPUTATIONAL PROCEDURES 

As in the steepest-ascent method, this method requires the use of a’ 
digital computer* The sequence of operations is summarized here* 

1. Select a nominal control history and initial values of v. This 
can be done by starting with a control program and V generated by the 
Steepest-ascent method, where 
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V 




from Eq. (35) • 

2. Integrate the system differential equations as in the steepest - 
ascent method. 

3. Integrate Eqs. (58), (59) > and (60), backward with the appropriate 
terminal conditions. If the determinant of M changes sign or H becomes 
positive, store the current value of the time and stop the integration. 

The reason for this is explained below. 

4. From Eq. (73) , compute the new control history while integrating 
the system equations forward. 

5- Compare the value of cp and >|r obtained to those obtained on the 
nominal trajectory. If the pay-off, cp, or the terminal constraints, 
have become worse the run is considered unsuccessful and a tighter bound is 
placed on 5u by increasing the magnitude of K. If , on the other hand, cp 
and ^ are the same or have improved, the run was successful and (3) and (4) 
are repeated. 

6. The program is terminated when no change occurs in the pay-off 

or the constraints and | K | « |h^| • At this point the control history is 
stored on punched cards or tape. 

7 . If in step (3) the determinant of M changed sign (this condition 
is called a "conjugate point" in the Calculus of Variations), or H 
became positive (which indicates the Legengre condition is not satisfied), 
the integration in (4) is begun at a slightly later time than this con- 
dition occurred. Normally, on subsequent iterations, this point will 
move backward to the beginning of the trajectory and disappear. 
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V. PROGRAM EXAMPLES 



A single-stage rocket trajectory problem as described below was 
programmed utilizing each of the methods discussed. The ALGOL computer 
language was used and the programs were run on a Burroughs B-5500 digital 
computer . 

Assuming the rocket is launched from an airless, non-rotating Earth, 
the state equations are 

= r = V sin (y) 

Xg = 6 = ^ COS (y) 



= V = g, 




o Iw, 



p sin (y) 

2 

r 



Xi = 




fj cos (y) , V COS (y) 



r 2 V 



where r = altitude measured from the center of the Earth , V = velocity, 
y = flight path angle, gQ = gravitational acceleration at the Earth T s 
surface, T = thrust (assumed constant), u = thrust angle (measured from 
the velocity vector), Wq = initial weight, Isp = specific impulse, t = time, 
\i = universal gravitational constant, 0 = downrange angle. The initial 
conditions are 

r(0) = R^ (Earth radius) 

0 ( 0 ) = 0 

V(o) = 100 ft/sec 
y(0) = 89*87 degrees 
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It is desired to place the payload in a circular orbit of maximum altitude* 



hence 



9 = x 1 (T) = r(T) 




Appendices (A) and (B) contain listings of the steepest-ascent program 
and second variation program respectively. Comments are inserted at strategic 
points which explain the sequence of operations. 

The steepest-ascent program contains logic for a single or dual stage 
rocket. It was rim in the single stage mode to generate a nominal control 
history for input to the second variation program and to compare results. 

It was also run in the tiro-stage mode to test the coast optimization logic. 

The input data for the steepest-ascent program are 

1. Initial velocity (f eet/second) (must be non-zero). 

2. Launch angle (degrees). 

3. Duration of the first stage burn (seconds)* for single stage 
rockets this quantity is the total burn time. 

b. First stage thrust (pounds). 

5- Second stage thrust (pounds)* for single stage rockets zero 
is input. 

6. First stage fuel flow rate (pounds/second) . 

7. Second stage fuel flow rate (pounds/second) * for single stage* 
any non- zero number. 

8. Rocket liftoff weight (pounds). 



2b 
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Second stage weight after separation (pounds) , for single 
stage j any non-zero number. 

10. Initial value of coast duration (seconds) , for single stage > 
zero. 

11. Duration of the second stage burn (seconds) , for single stage., 
zero . 

12. Coast weighting factor, A. 

13. Number of stages (l or 2). 

The input data for the second variation program are: 

1. Initial velocity (feet/second) . 

2. Launch angle (degrees). 

3. Duration of rocket burn (seconds). 

4. Thrust divided by initial weight (pound/pound) . 



7 • Integration step size and data storage interval. 

8. K (See Section IV-A) . 

9* Nominal control history. 

Other parameters which the user may desire to change must be changed in- 
side the program or incorporated into the READ statement. 

For the single-stage run s, the following input values were used: 




Launch angle 



Initial velocity 



100 ft/ sec 
89*87 degrees 



Duration of burn 



220 seconds 



Fuel flow rate 



Thrust 



430,000 pounds 
1433.3 pounds/second 
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Liftoff weight 



333,770 pounds 
678,914 
-93.58 



V 1 




V 


-93.58 


Step size 




K 





For the two- stage runs, data for the Gemini -Titan II system were used: 



Initial velocity. ...... 




Launch angle 




First stage turn 




First stage thrust 




Second stage thrust .... 




First stage fuel flow.. 





Second stage fuel flow 327 *7 pounds/second 



Liftoff weight 




Second stage weight.... 




Coast duration 




Second stage burn 




Coast weighting factor. 


0.1 


Number of stages . . 0 ... . 
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VI. RESULTS AND DISCUSSION 



A. STEEPEST-ASCENT VS. SECOND VARIATION 

When the steepest-ascent program was run in the single-stage mode, 
the nominal trajectory attained an altitude of 196,015 feet. The errors 
in meeting the terminal constraints on the velocity and flight path angle 
were 426 feet per second and 1.146 degrees. On the third iteration the 
altitude was improved to 260,427 feet with terminal errors of 0.67 feet 
per second and 0.025 degrees. The program was terminated when the desired 
altitude change , dcp, became less than 5^000 feet. At this point fifteen 
iterations had been completed. The terminal altitude achieved was 318,126 
feet with terminal errors of 0.64 feet per second and 0.003 degrees. The 
associated control history was punched on cards and values of and 
were printed out. The program ran five minutes and nine seconds. It is 
estimated that this time would be about halved if the program were run on 
an IBM 7090 computer. 

The output generated in the steepest-ascent program was used as input 
to the second variation program. As expected, the nominal trajectory 
attained an altitude 318,126 feet. On the succeeding forward integration 
the trajectory was totally unreasonable. The control deviations were 
made smaller by increasing the initial magnitude of K but this failed to 
improve the results. Small variations were made in the input values of. V 
which caused relatively large changes in the results but it was not clear 
how to make adjustments which would improve the performance of the program. 
The rimming time was far in excess of the steepest-ascent program, taking 
over three minutes to compile and complete just one iteration. 
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B. TWO STAGE ROCKET TRAJECTORY WITH COAST OPTIMIZATION 



As stated in Section IV, the input data for this problem were those of 
the Gemini-Titan II launch system with an arbitrary choice of ten seconds 
for the initial coast duration. Lambda., the coast weighting factor, was 
chosen as 0.1 since earlier runs indicated that a value of 1.0 caused the 
coast variations to be insignificant. This choice proved to be satis- 
factory. 

The nominal trajectory attained an altitude of 392,564 feet with 
terminal errors of 25 feet per second in velocity and 0.84 degrees in flight 
path angle. On the first iteration, where the pro gran attempts only to 
meet terminal constraints, the altitude was improved by l8,Ol8 feet, the 
terminal errors were 0.09 feet per second and 0.0025 degrees. On the 
fourth iteration the coast time was reduced to eight seconds. At this 
point, a new nominal for the portion of the trajectory following first 
stage burnout was computed using the new coast period. The result of 
this change in coast was an improvement in the velocity constraint of 
12 feet per second, and a degradation of the flight path angle by 0.25 degrees. 
The terminal altitude achieved on this iteration was 573,450 feet. On the 
fifth iteration the coast period was reduced to two seconds, this accounted 
for an altitude improvement of 30,973 feet. On the sixth iteration the 
coast period was reduced to zero, this improved the terminal altitude by 
13,477 feet. In both instances cited above where the coast period was 
changed, the terminal constraint errors were diminished. 

The altitude improvement specified on the second iteration was 50^000 
feet. The program was allowed to run until this figure was reduced to less 
than 1000 feet. This proved to be rather wasteful as the terminal altitude 
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failed to improve significantly after the desired improvement was reduced 
to 6,250 feet, which occurred on the eleventh iteration. In Fig. 1 the 
terminal altitude, corrected for terminal constraint errors, achieved on 
each iteration is illustrated. It clearly shows that twelve iterations 
were sufficient to converge to the optimum. 

Figure 2 illustrates the convergence of the control history. The 
discontinuities which occurred are due to using discrete feedback at 
twenty second intervals rather than continuous feedback. The curves are 
of different lengths due to the changes in coast period which occurred. 

Figure 3 illustrates the action of the terminal error control scheme. 
Each time the errors became excessive they were reduced to essentially 
zero in one iteration. 

The last terminal altitude achieved was 6l6,573, a 57*2 per-cent 
improvement over the initial nominal. 

C. SENSITIVITY OF THE STEEPEST -A SCENT PROGRAM TO CHANGES IN INITIAL 
CONDITIONS 

The sensitivity of the steepest-ascent program was tested to determine 
the capabilities of the convergence scheme. Lauch angles of 89.87 and 89*95 
degrees were tested on the single stage trajectory. The resulting terminal 
altitudes were 196,015 and 803,^28 feet. In the latter case the terminal 
errors were 1858 feet per second in velocity, and bO.b degrees in flight 
path angle o In both cases the program converged in twelve iterations to 
approximately the same optimum altitude. 

An attempt was made to solve the problem given a ninety degree 
launch angle which failed. Further testing at % lower launch angles was not 

accomplished due to computer time limitations, however, it is believed that 
the tests conducted amply illustrate the virtue of the method. 
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VII . CONCLUSIONS 



The terminal control feedback scheme due to Bryson and Benham and 

rf B 

the method due to Rosenbaunr of leaving the pay-off unconstrained in order 
to satisfy terminal constraints have been combined and altered as necessary 
to produce a simplified steepest-ascent optimization procedure. This 
procedure has been shown to be successful in solving a typical rocket 
trajectory problem including optimization of an interstage coast period. 

In this closed- loop procedure the change in control is computed at 
each point utilizing continuous or discrete information on the deviation 
from the previous trajectory. It is closed loop because the program 
continuously (or periodically) checks on how it is doing in its attempt 
to satisfy the terminal constraints. The advantage of this procedure is 
that larger deviations from the nominal can be tolerated while still 
maintaining a bound on the terminal errors. It is, therefore, possible 

to move rapidly toward the optimum trajectory as is illustrated in Fig. 2. 

6 

The second variation method due to Bullock has been applied to the 
same problem with questionable results. Consultation between the author 
and Mr. Bullock has failed to uncover possible flaws in the theory or the 
programming technique. Bullock has shown in several examples that the 
method is successful, however, in each example the Hamiltonian did not 
depend explicitly on time. In the rocket problem this is not the case. 
Although this fact has no theoretical bearing on the problem, it is the 
only major difference between Bullock's examples and this problem. 

Experimentation revealed that the rocket problem is extremely sen- 
sitive to the choice of v. This was not the base in the examples 
presented by Bullock. It was thought that the values of V generated by 
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the steepest-ascent program would be very close to correct , however , in 
view of the failure , this premise was laid open to question. Hence, one 
reasonable conclusion that can be drawn is that v was chosen incorrectly 
and that there is no intuitive or analytical method presently available 
to make the proper choice. There exists, of course, the possibility of 
error, but considerable time and painstaking effort has been expended to 
minimize this possibility. 

Failure of the second-variation method notwithstanding, some con- 
clusions may be drawn with respect to the advantages and disadvantages 
of the two methods . 

1. Understanding the theory involved in the second-variation method 
requires a background in the Calculus of Variations , whereas the steepest- 
ascent method presented does not. 

2. The second variation method apparently requires a good estimate 
of the sensitivities , v ; while the steepest-ascent method only requires a 
guess of the nominal control and will tolerate a fairly poor guess. 

3. The second variation method requires backward integration of the 
2n2 + 2n equations M, N, b, and A (where n is the number of state variables) 
and a matrix inversion at every step of the forward integration. The 
steepest-ascent method requires backward integration of $ and J which is 
less than 2n^ equations since J is symmetric. The matrix inversion can 

be done at less frequent intervals using the sampled data feedback method 
suggested. The second variation method thus requires more computer time 
per iteration. 
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VIII. ADAPTATION OF THE PROGRAMS TO OTHER TYPE PROBLEMS 



An attempt was made to generalize each program so that they could he 
easily adapted to other problems. This required a large number of sub- 
scripted variables and matrix multiplication loops in the subprograms 
containing the differential equations. Since the subprograms are called 
twice per integration step they must be as efficient as possible. The use 
of subscripted variables and loops is most inefficient and results in more 
than doubling the running time. 

Adaptation to another problem is still relatively simple however. 

The programs contain sufficient comments to indicate where necessary 
changes must be made. 
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APPENDIX A. STEEPEST -ASCENT COMPUTER PROGRAM 



RE.GIN COMMIT ERNEST C. LUDERS PtOX 219 S T E E PE S T - A SC E N T J 
INTEGER FLAG#FA1L#CUUNT#QUIT> 

REAL T1#C0m$T#FF 1#FF2#ISP1#ISP2>WR1#WR2#W1#W2»T2#MM1#MM2# 
RE/ MASQAT, lUVERWO# GC# K# ISP# TF# HH# SAMPLETIME# 
VO#PREOhFG#PREHFO#CHFC#HFO#Jl2#J13#PREVF#PREGAMF# 
PrttHF# I.) V F # DGAt.F# ROUND# BCUN'02# UET# EX# CNR# P# 
Mm2,LAMOA»SUM#ULDCOAST#TB2#G#DHF# 

INTEGER MAXir;DEX#L#M#I#lTER#STAGES#OLDHAXlNDEX# 
COMPUTECOASl/COASrCOHPUTE; 

REAL ARRAY D [ 1 : 3# 1 : 3 3 # A [ 1 : 3 3 * Y P C 0 : A #0:1103# 

SAVE A R n A Y xP # OLDXP I 0 : A # 0 ’• 200 3 # LP C 0 ! 22# 0 : 200 3 # ERR[0:33# 
TEMPI 0:43, D E L C C : 3 3 # C HK [ l : 3 # 1 : 3 1 , 

P-'PU :2# 1 S23# PHI #OLCPHI [ 0 :200 )#uEI# 
j 1 n V [ 1 : 3 # 1 : 3 3 J 

LA3EL Li# 12# L3# LA, L5 # L6 # HELL# 

FORMAT f- 1 (//# X 25 # "AT TIME T = "# 13# " THE OETERMl", 

” NANI OF J = "# 

E18.11# //# X 3 3 # 3E23.U# //# X13# "CHECK MATRIX = ”# 

3E23.11,//# v 3 3 # 3E23.11# //)# 

F 2 ( X 2 0 # "THIS RUN REQUESTED AN ADDITIONAL:"# //)# 

F 3 (X2?#F12.2# ” FEET OF ALT ITUde" # //# X30# F9.2, 

" FEET/SEC UF VELOCITY"# //# X32# Fll.6# 
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" degrees nr flight path angle"# //)# 

FA (X20# "THE RUN ACTUALLY YIELDED AN AOC I T I ORAL : "# //)# 
F5 ( X 2 0 # "TERMINAL ALTITUDE = "# F10.2> 

" FEET"/ / # X 2 0 # "TERMINAL VELCCITY"#X17#"="# 

X 1 # F 9 . 2 # " FEET/SEC"# /# X20 # 

"TERMINAL FLIGHT PATH ANGLE = "# F11.6# 

" DEGREES"# / # 

X20# "ORBITAL VELOCITY AT THIS ALTITUDE = "# XI# F9.2> 

" FEET/SEC")# 

F 6 (X2C#"C0AST DUKA 1 ION"# X2C# " = "#F 10 , 2# " SECONDS")# 

F 7 C X20# "CORRECTED TERMINAL ALTITUDE ="#F10.2#X7# 

"FEET")# 

HIS1uRY(4F13,12)# 

DEFINE LOOPL = FUR LH#2»3#A DO X# 

FILE GUT CAriQS 0 (2#10)# 

COMMENT MATRIX MULTIPLICATION PROCEDURE GOES HERE# 



COMMENT MATRIX INVERSION PROCEDURE GOES HE RE x 



COMMENT THIS PROCEDURE CONTAINS THE SYSTEM DIFFERENTIAL 
EQUATIONS Ai\Q IS USED ON THE FIRST NOMINAL# ON RUNS IN 
WHICH THE COAST THE HAS BEEN CHANGED # AND QN THE FINAL 
PRECISION HUN# 
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PROCEDURE FU\CTCTj> X , F ) ; 

value t; real t; real array x, ecu; 

BEGIN REAL H/ V > S, C^FEE>TS; INTEGER i* j; 

COYMENT IinTlRPOlATE in I HE CONTROL array for the 
PROPER valul; 

I <- r-N r IER c T/HH) ; J <- IF I<MAXINDEX then I + l ELSE I 
FEE ♦ t«3<-PHI[I]) + (T/HH - I)x(PHI[J] - Q ) > 

r <- x[ 1 1 + re; 

s <- s i n ( x ( a i ) ; c «• c o s ( x t a ] ) ; v «- x c 3 1 ; 

TS«-T-T2; COMMENT TIME FROM STAGE TWO IGNITION; 

IF T < T 1 OR STAGES =1 THEN BEGIN 
TGV£RrtO<-FF t / rt 1 ; 

MASRA1 «■ 1 - TOVERWOxT/ ISPi; 

ENO ELSE tuverwo<-c; 

IF T > T 2 AND T2/TF THEN BEGIN 
I 0 V £ R « 0 <- F F 2 / w 2 ; 

MASRAK-1-T0VERW0XTS/ISP2; ENO; 

COMMENT THESE ARE THE SYSTEM DIFFERENTIAL EQUATION'S; 

F [ 1 I <• Vxs; 

F [ 2 I <- VxC/R; 

F [ 3 ] «• (Q <• GOXTUVERWO/MASRAT JxCOSCFEE) - Kx$/R*2; 

F [ 4 J f -K/Vx C/R*2 + QxSIN(FEE)/V + F [ 2 1 ; 

end funct; 



COMMENT THIS IS the BACKWARD INTEGRATION PROCEDURE WHICH 
CONTAINS THE ADJOINT EQUATIONS AND THE CONTROLLABILITY 
MATRIX EQUATIONS; 
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PROCEDURE LINBAK (Td> LS* LF); 



value tb; 



real r&; real array ls, lfcii; 

begin REAL R> V, S> C> T# SF r l.l> L2, L3> L A > IM5 
INTEGER I>J; 

real ts; 

REAL Lb»L6>L7>L8>Ni,N2*N3,M>N5>N6>N7»N8; 

REAL TSM,TCM>M1,M2,M3>FEE>L13>LH#L33#L34>L43>L44; 
COMMENT T IS BACKWARD RUNNING TIME; 

T «- IF - TO; 1 «■ ENTIERCT/HH); 
if iso then i <- u; 

J <- IF I<MAXINDEX THEN 1 + 1 ELSE I; I N T «- T/HH - U 
Ll<-LS[li; L2«-LSC2]J L3«-LS[3i; L'4«-L$U]; 

L5*-L S [ 5 ] ; L6*LSl6]i L7<-,LS[7j; L8<-LSC8]; 

M <- L S [ 9 J ; *42«-LStlOi; N3«-L Still* N4«-LS[l2i; 
■N5«-Lo[l3i; N6«-Lall/4i; N7«-LStl5]; N8*LSCl6i; 

CUN'MENT INTERPOLATE FUR THE PROPER VALUES OF THE STATES; 

R <- (Q<-XP C1>11) «• INTxtxPt \> J3 “ 0) + RE; 

V <- (Q«-XPC3#IJ) + INTx(XP[3*J] - Q ) ; 

S *• SlN((0<-XPtA,I]) * INTx(XP[4j>J] - Q ) ) ; 

C <• COS( ( Q«-XP I <\ f 1 ] ) + INTxCXPCA/J] - Q ) > ; 

SF <■ S I N ( £ Q* PH I C i ] ) + I N T x ( P H I [ J J - 0)); 

TS<-T-T2; CUKMENT TIME TJ GO UNTIL STAGE TWO IGNITION; 

IF T > T 2 AND T 2 7 IF THEN BEGIN 
T O.V E R ^ 0 «- F F 2 / W ? * 

KASRATM-TUVEHWOXIS/ISP?; end 
ELSE TuVERW 0 «- 0 ; 

IF T < T 1 OR STAGES =1 THEN BEGIN 
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tdverwo«-ffi/wi; 
MASRAT*i-TOVERw ox t / i spi; 



ENOJ 

COMMENT T H E 3 E ARE THE ADJOINT DIFFERENTIAL EG NATIONS; 



IF [ 1 ] *■ 


•VxCxL?/R*2 +2xKx$xL3/R+3 - Cx(V - 2 x K / ( V x R ) ) xL a / 
R * 2 ; 


LF [ 2 ] «■ 


o; 


LFt3] «■ 


+ Sxi_l +C*L2/R + C x ( l + K/C V*2xR ) )xL4 /R - GO* 
T0VERw0xSF/(V*2xHASRAT)xLAJ 


LF[4] <• 


+ V x q x l 1 - V*S*L2/R -KxCxL3/R*2 - Sx(V - K / ( * R ) ) 
xL4/R; 


L F [ 5 ] «■ • 


'V x CxL6/R*2 + 2x ><xS x L7/R* 3 - Cx(V - 2 X «/ ( VxR ) ) xL3/ 
R*2T 


LFC6] «■ 


o; 


I.FI7 J <• 


+ 3 * L 5 + C x L 6 / R tcx(l + K / C V * 2 x R ) ) x l 5 / R - GO* 
T 0 V E R N 0 x S F / C \/ * 2 x M A S R A T ) x |. 8 ; 


LF c e ) f 


+VXC X L5 - V X S X L6/R -KxCxL7/R*2 - Sx(V - K/(VxR)) 

xl8/r; 


LF [ 9 ] «■ ' 


■V x Cx-M2/R*2 + 2xKx$X|\3/R*3 - Cx(V - 2 *K/(VxR))xNA/ 

r*2; 


L K C 10 3 <• 


o; 


LFC11] «■ 


fS*Wl + C x i J 2 / R + C x ( l + K/(V*2xR) )xN/j/R - GO* 
r 0 V t‘ R H 0 X 3 F / (. V * 2 X M A s R A T ) X N '4 ; 


L F t 1 2 ] <■ 


+ V x C x N 1 - V*sxN2/R ”K*CxN3/R>2 - Sx ( V - K / ( VxR ) ) 

•X N A / R ; 


L F C 13 3 «■ ' 


•\/xCxN6/R*2 +2*KxSxn7/R-*3 f C x(V - 2 x K / ( V x R ) ) x N8 / 

r*2; 


LF C 1 A 3 «- 


0 * 



IF [ 1 5 3 «■ +SxrJ5 +CxNo/K + Cx(1 + < / ( V * 2 xR ) ) x N 3/ R - GO* 
T0VertW0xSF/(V*2xMASRAT)xN3; 

L F [ 1 6 ] <- + V x c '< N b - VXSXN6/K -KxCxN7/R*2 “ Sx(V - K/CVxR)) 

< N 8 / R / 

FEE «• (Q^PHltU) + IUIX(PHICJ]-Q); 

TSM *■ GOxTnVER»iO/HASRATxsrKKF.E); 

T C • 1 «• GOxTlWERWO/MASRATxCQSCFEEJ/VI 
Ml «• " L 3 x T 5 M + 1 4 * r C H > 

M2 *■ -N3XTSM f N a x r C M } 

M3 «■ - N I x T S M + N 6 x r C M i 

COMMENT THESE are the controllability matrix EQUATIONS; 
LF'ClFJ <- Ml*?; L F I 1 8 3 «■ M2*2; LFC193<- M3*2; 

L F [ 2 0 J <- M1XM2J L F [ 2 1 3 «- M 2 x M 3 : LF[ 22 3«- fOx^i; 

end linbak; 



COMMENT r H I S IS THE FORWARD INTEGRATION PROCEDURE WHICH 

COMPUTES the new control; 

procedure slope c t, x, ft; value t; 

REAL T; REAL ARRAY X, F [ l 3 > 

BEGIN real R , V, S, C, T3, I ND , BK I NO/ FEE , TSM, T C M , Ml, 
M2, M3, TS; 

LABEL LI, L2, L3> LA, TCQCL0SE»Lb,L6» INTEGER I> j; 
TS<-T"T2; CUMMENT TIME FROM STAGE TWO IGNITION; 

IF T < T 1 OR STAGES =1 THEN BEGIN 
T 0 V E R w 0 4- F r i / w i ; 

MASRAT*-l-TQVEHWOxT/ISPi; I$P«-ISP 1J 
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END ELSE TUVtRWO«-0; 

IF T > T 2 AMU T 2 T F THEN BEGIN 

T0VE^rt0«-FF2/N2; ISP<-ISP2/ 

masrat<-i-tuveufuxtb/isp2j end; 

IF T >yr)UNOx • 999999999 UR COMPU TE C 0 A S T = 1 THEN 
begin 

COMMENT THIS IS THE SAMPLED DATA FEEDBACK SECTION AND IS 
HIT AT THE SAMPLING INTERNAL AM) AT T2 WHEN A NEW COAST 

period is computed; 

IF COMPUTECUAsTM THEN 
begin 

I N 0 «- 6 nuND/HH; 

HOUNDS BOUND + 3 AMPLE T I MET 
End else IND<-ENTIER( T2'/HH+. ooooooi ); 

B r\ I NO <- M A X I N BE X - I N D ; SUM<-0; 

COMMENT SINCE J GETS SINGULAR NEAR TF » SKIP THIS SECTION 
IF T iS CLOSE Tj TF; 

IF BOUND > TF-30 THEN GO TO TOOCLOSE* 

FOR I <• 1 > 2 > 3 00 JE I C I » I 3 «■ JIHVCI# I )<-LP[I + 16/0KlNO) + 
(IF T > f 2 x , 9 9 9 V 9 9 OR COMPUTECOA ST = 1 THEN 0 [ I > I 3 / 

L A m 0 A ELSE 0 ) ; 

J c. I [ 1 >21 <- JIN VI 1>23<-JEI 12 / 1 3«-JINV[?> 1 ]«-LPt20-»BKlNo3 + 
(if T > T ? x , 9 9 9 9 9 9 OR C OHPU T E C 0 A ST = 1 THEN 011*23/ 

l a m o a else o); 

JET(2*3]<-JiNVC2*3]t-JEIC3>2]«-JINV[3*2]<-LPC2l/0KIND] + 
(IF T > T 2 x , 999999 OR COMPU TECO AST* l THEN D[2*33/ 

L A m 0 A ELSE 0 )p 

UtIU,33«-JlNV[l,33<-JEIt3>13<-JINV(3,l]<-LP[22>BKlN03 + 

b2 



(IF T>T2*. 999999 OR C OMPUTECQ A ST= 1 THEN DC1/3J/ 
LmMOA ELSE 0); 

IF C 0 M P U T E C 0 A S T = 1 THEN C G M P U T EC U A S T «• 0 / 
p <- timecd; 

IMERTC JINV, 3/ I); 

If- 1 = 1 THEN BEGIN 

WHITECO'THE j matrix is singular at tire = "/ 

F6.2>/ T); 

CURRENT IF THE j MATRIX BECOMES SINGULAR TERMINATE THE RUN/ 
GU TO HELL END; 

LUOPL TEMPtLI «■ OLOXPtL/ INQ1J 

LUOPL SUM <- SUM + LF [L/ 0 K I NO I x ( X C L 3 - TEMPtLI); 

oelcu *■ dht - sum; sum <■ o; 

LUOPL SUM <• SUM+LPtL+S/ 8 K I NO ] x ( X C L 3 - TEMPtLI)/ 

ut.Lt 2 J *■ dvf - sum; sum <• o; 

LUOPL SUM «■ SUM + L p [ L + 12 / BKIND)x(XtL) “ TEMPtLI); 
DELI 31 *■ DGAMF •* SUM; SUM *■ 0 } 

EnRtn * err 1 2 ] <• err 1 3 1 «■ o; 

FUR L «- 1 / 2 / 3 DU FUR DO 

EKRtL] «■ F.RRtLI + JlNVfL/ M]xOEL[Mj; 

CURRENT ERR IS THE MATRIX OF LAGRANGE MULTIPLIERS/ THIS 
IS THE END UF THE SAMPLED DATA FEEDBACK SECTION; 

TUQCLOSE: ENrj; 

COMMENT THE CONTROL IS COMPUTED BETWEEN HERE AND LA/ NOTE 
THAT THE CUnTROL IS COMPUTED ONE STEP AHEAD OF THE CURRENT 
TIRE TO PERMIT INTERPOLATION/ 

LIS IF T £ 0 0 IJ N D 2 x , 9 9 9 9 9 9 9 9 9 THEN 



BEGIN IF T = 0 THEN BEGIN I<-o; 



GO TO L2 END/ 



L3: 1 <- 90UND2/HH + 1) B0UNn2 * 80UND2 + H H ) 

IF iJ0UNi)2>TF I HEN GO Til 1. 4 ) 

L 2 : FEE «■ PHICIJ) V <- X F [ 3 , I ] ) BKINQ «■ MAXlNQEX - I 

TSM30UND2) 

IF 3 0 U 'N* D 2 > T 2 I HEN T5t-B0UN02-T2; 

TSm«-GO x r IJVCPNO/ (0<-( 1-TOVEPW OXTS/ ISP) )XS INC FEE)) 
ICrt «■ GOxTUVEKWO/QxCCSCFEE )/V) 

Ml <- -TSMxlP L3^ ; 3KINO) + TCMxLP [4/1KIN0]) 

M 2 <- -TSMxlPC 11 ,'TKINQ I + T C Mx L P [ 1 2 > B K I NO ) ) 

M3 «• - TSMXLPC 15V3KIN0 ] + T C Mx LP C 1 5 , 3 K I N D I ) 

COMMENT CUMHUTE THE NEW VALUE GF THE CONTROL) 

PHI C I ] «■ (JLOPHI C I ] «- F H I [ I 1 ) +H1XERRC1) +M?xERR[ 2] 
+ M 3 x E R R [ 3 ] ) 

Ir UEhTIERC T 2/HH+ .OOOOOOl )+l AND STAGES=2 THEN 

CUMPUTECUASTM 

ELSE COMPUTCCUASff-O) 

IF T = 0 AND 1=0 THEN GO TO L3 END) 

L A $ I «• ENT I ER ( T/HH ) ) J <■ [ F I<MAXINDEX THEN 1 + 1 ELSE I 
FEE <■ CCX-PHICD) + (T/HH - I)x(PHI[J] - Q)) 

R <• a [ 1 ] + RE) 

S <- S I N ( X C A I ) ) C <- C 0 S ( X [ 4 1 ) ) V «• X [ 3 I ) 

COMMENT HERE ARE THE SYSTEM DIFFERENTIA!. EQUATIONS) 

F[ll <- V •< S ) 

F [ 2 I «• V < C / R ; 

F [ 3 ] <- CO «■ G 0 x T U V E R R 0 / m A S R A T ) x C 0 S C F E E ) - K*S/R*2) 

F [ 4 J <■ - 41 /Vx C/R*2 + QxSlN(FEE)/V + FC2J) 

COMMENT IN THIS SECTION THE NEW COAST DURATION IS COMPUTED) 
IF CUMPUTECUAST=1 UR C 0 A S TC 0 MP U T E = 1 THEN GO TO L5 



else go ro L6; 

L5: COAblCOMPUTL<-CJASTCOy,PUTE + i; 

A [ 1 X-Q* 

A[2]<--Kx$/R*2 * (-[3]; 

A[3 I*-K/V*C/R*2 + F [ 2 ] - F C 4 I * 

BKINOf MAXI NOEX-EN TIER(T2/HH+. OOOOOOl ) } 
MM1«-A[2]xLP[3/3K1ND] ♦ A C 3 ] x L P [ A * 8 K I NO ] ; 

M M 2 <• A [ 2 I x L P l 1 1 > 3 K I .\| 0 ] + A [ 3 ] xL P t 1 2 * 0K I NO I f 
MM3«-A[2 ]xLPC 15, HK INO] + A [ 3 1 xL P [ 1 6 , 8K I NO ] ; 

D [ 1 * 1 ] <• MM 1 * 2# 0 [ 2 > 2 ] «■ MM2*2i 0 [ 3 > 3 D <• MM3+2) 

D [ 1 > 2 J <- D [ 2 , 13 <- MMIXMM2; 

D [ 1 * 3 ] «• 0 C 3 j. 1 3 *■ M -1 1 x M M 3 ,* 

D [ 2 * 3 ] <■ 0 [ 3 » 2 ] <• M M 2 x M s' 3 ; 

IF CUMPUTECOASTal THEN GO, TO L 6 * 

COAST <- ( 0 L 0 C 0 A S 1 «- C 0 A S T ) + (Q<-(VMl x ERR[ 1 ]+MV2*ERR[2 ] 
+MM3 *£Rh[ 3 J ) /LAMi)A)J 

COMMENT COAST is set TO the nearest integer divisible OY 

THE STORAGE INTERVAL* THIS AI0$ IN KEEPING INDICES STRAIGHT 
i\R I It ( <"0F.L T ACUAST = "*£20. 1C>,Q); 
C0ASWE\ITI£R((C0AST + .9)/2)x2; 

if cuast<o or stages =1 then coast<-o; 

WRITE(<" COAST = ",F2C.10>* COAST)* 

Lo: END SLOPE* 

COP PENT INTEGRATION PROCEDURE GOES HERE; 

COMMENT MAIN PROGRAM BEGINS HERE* 

COMMENT INITIALIZE BELLS, FLAGS* AND CONSTANTS) 
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ITEft<-0; FLAG<-0/ F A I L <- 1 : COUNT *-0/ GUI T*o; 



RE «• 2.09«?7; GO <- 32.17; K <- GCxRE*2) 

COMMENT SAMPLETIME IS I he feedback sampling interval and 
HH Is THE data STORAGE INTERVAL* HH IS ALSO THE INTEGRATION 

STEP size; 

sakpletime «■ 2C ; hh <- 2 ; 

Lis REA0(XP13/C]/XP[4/0J/T1/FF1/FF2/NR1/WR2/W1/W2/CCAST/TB2/ 
LAK0A#STA5ES)[L4]; 

COMMENT THE input DATA IS 





XP [ 3 / 0 ] = INITIAL VELOCITY 
X P [ A / 0 ] = LAUNCH ANGLE 
T 1 = FIRST STAGE BURN TIME/ 


COMMENT 


F F 1 = FIRST STAGE THRUST 

F F 2 * SECOND STAGE -THRUST 

rtRi = FIRST STAGE FUEL FLO K RATE/ 


COMMENT 


w R 2 = SECOND STAGE FUEL FLOW RATE 
W1 = LAUNCH WEIGHT 

W2 = SECOND STAGE WEIGHT AFTER SEPARATION/ 


COMMENT 


COAST = INITIAL CHOICE CF THE COAST DURATION 
T32= SECOND STAGE BURN TIME 
LA MG A = COAST WEIGHTING FACTOR/ 


comment 


STAGES = NUMBER OF STAGES; 



0LDCQAS1 vCOAST J 
T2<-Tl+CLAST; 

T F <■ T 2 + T b 2 / 

0 L D M A X I N D E X f Mi A X I N D E X <- E N I I E R ( T F/ HH+ . 5 1 ) 

ISPUFFl/hRlJ I SP2<-FF2/WR2; 

XPtl/OJ <- X P [ 2^03 «■ o; XP[4/0] <- XP[ A/03/57 ,2957795 1 3 1 ; 
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COMMENT GENERATE Or. READ IN THE NOMINAL CONTROL HISTORY; 

FOR J<-G STEP 1 UNTIL MAXINDEX DO PHICII <■ 0> 

COMMENT INITIAL CONDITIONS FOR BACKWARD INTEGRATION; 

FOR J <- 2 STEP 1 UNTIL 22 00 LP[I>C]<-OT 

FOR I>1 STEP 5 UNTIL 16 QU LP[I#C]M> 

comment compute the nominal trajectory; 

A D A M S ( A » hH, T F ^ HH > T F , 0/ C , XP, FUNCT); 

COMMENT perform thf. backward integration; 

L 2 ; ADAKS(22,HH, 0/ TF, HH* TF » 0/ 0> LP > LINBAK); 

COMMENT store THE OLD VALUES CF THE STATES BEFORE COMPUTING 

a nek forward trajectory; 

FOR I <■ C STEP 1 UNTIL MAxiNCEX DO LOOPL OLDxF[L>I] *■ XP[L>II 
COMMENT BETWEEN HErE AND L 3 THE ALTITUDE CHANGE DUE TO 
TERMINAL ERRORS IS COMPUTED; 

PREHF «■ XPE WMAXlNDEXi; PREVF <■ XPt 3> MA X I NDE X ] ; 
pregamf «- xpi AjDaxIndex]; 

J12^LP[20/MAX1NDEX]; 

J13«-LP[Z2»MAXlN0EXi; 

TM'Pt i/ ik-lpl lb^ max index]; 

!MPt2/i]«-TMP[l»2]«-LPt21/MAXlNDEX]; 

1MPL2*2] < *LP[ IV / MAXINOEX]; 

I N V E H T ( T M P > 2 r I) ; 

IF 1 = 1 THEN BEGIN WKlTE(<"THE LAMBDA MATRIX IS S I NGUL AR"> ) ; 

GU TO LA END! 

lHP[t^l]«-Jl?xlMP[l,n+J13xTMPC2#U; 

Tf'PC 1>2 3<-J1?X |MP[ 1,2 34 J13XTMPC2/21; 

VU <- SCRKK/CPREHF + RE)); 

PRE0HF04-TllP[W13x(VU-PREVF)-TMPCl>23xPREGAMF3 
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P*EhFG<-PHt HF +PREDPFQ; 

rtKIUC<”OLD CORRECTED HF = " » F 2 0 . 1 0 > , P R EHF 0 ) > 
L 3 t IF FLAG =0 THEN ELGIN 

dhf<-frldhfo; 

ovi<-sgr7(k/(prehf + uhf + re))-prevf; 
dgamf< -pregamf; end; 

BOUND «■ ROUND 2 <- COASTCOMPUTE <- COMPUTECOAST ♦- o; 
COMMENT COMPUTE THE NE*v CONTROL WHILE INTEGRATING THE 
SYSTEM EQUAIIONS FORWARD; 

ADAPSCA, HH, 0, T F , HH> IF, 0, 0/ XP, SLOPE)! 

COPMENT THE FOLLOWING LOGIC CAUSES THE CONTROL PISTCRY 
TO BE PRINTED OUT EVERY FOURTH ITERATION; 

IF ENTIER(ITFR//!)xa = ITER THEN BEGIN 

.(RIIECFOR 1*0 SILP 1 UNTIL HAXINDEX DC PHICI3); 
WRITE! [PAGE! ) END; 

ITER*I TLR + i; 

0 L 0 h A X i N D E X <■ m A X ] N D E X ; 

IF COAS UULDC0AS1 ThEN 
BEGIN 

comment adjust storage locations OF PHI AND Fly new 

NOPINAL TO PROVIDE PROPER INPUT OF STATES FOR BACKWARD 
INTEGRA I IUN» 



.T 2 «- II + C 0 A S T ; 

TF*-T2 + TB2F 

0 L 0 M A X I N D E X «• f , A X I N D E X ; 
MAXINDEX<-LNTJER( TF/HH+.5J ); 
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L«-ENlIEN((COAST“ULDCOAST)/HH+.OCOOOOnJ 

F 0 K 1 <-LNTIER( ( Tl + OLUCOAST )/HH+ .0000001 ) STEP 1 

UNTIL OLDMAXINDEx DO BEGIN 
PulU+LX-PHltl)! 

U I<[M I EK( TP/HH+ . 000000 1 ) "1 THEN PHIU + 1]«-0J 

end; 

i«-eniilr(Ti/hh + .0 00000D; 

LOOPL YPCL/0]«-XPCL»I]! 

a D A f 1 S ( 4 7 HHr Tl> T F > HH/ T F/ 0/ 0/ YP> F U N C T ) ! 

FQR I-.4-I $TL P 1 UNTIL MAaINDEX DO 
LOOPL xP[Lf M] ayPIL>M-I ]; 

eno; 

COMMENT PRINT OUT THE RLSUlTS OF THE ITERATION! 

WRITECF2)! w H I TECF 3 > DHF/ I3VF* DGAMFx 5 7.2957/95131 )1 
NRITE(F-J); wRITLCFi/ Xp [ 1 > «AX INDEX ] - PREHF, 

XPC 3/HAXI nOE a ] - P R E '/ F > ( XP t A , MA X I ND£ X ] - PREGAM F) 
*57-295779513)! 

WRITE l f 5 * XP[1, MAX INDEX]/ XP[ 3/ M AXINDEX]/ 

APC«/MaXIN0E*Jx5/.2957/95131> SQRT(K/(RE + 

XPC 1 / MAX I NDEX J ) ) ) ! 

WRITECI- 6 , COAST); 

COMMENT COMPUTE THE CORRECTED TERMINAL ALTITUDE! 

V J «- S 0 R T C K / C X P C 1 / M A X J NL) E X ] t R E ) ) ! 

U n F 0 <• T M P l 1/ 1 )*( VO-XPC 3/MAXINDEX) )- 
TMPC 1 / 2 ] XXPC A / MAX I NO EX ] ! 

H F 0 <■ X P C 1 / :T A X I N 0 E X J + DHFO! . 

COMMENT PRINT OUT THE CORRECTED TERMINAL ALTITUDE! 

’’Y h i t e ( r 7 7 h f o ) ; 
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IF 0 *j I T = 1 THEN GO TO L6* 

COMMENT TEST TERMINAL ERRORS WITHIN TOLERANCE* 

IF AOSlOVF«-CVO-XP[i^MAXINOEXI))> t jO 

OR ABSCDGAmF«-(XP[ 4#MAXINOEX3)) > .02 
OR SIGNCDHF )/SlGN(XPCl,MAXIXOEXl-PREHF) then 

bEGIN 

fail«-faii. + i; 

C u U N T «- C 0 U N T ♦ 1 ; 

comment test for improvement in satisfying the terminal 

CONSTRAINTS OR IMPROVEMENT IN THE CORRECTED TERMINAL ALTITUDE* 
IF (A3SCDVF)<AdSCS«RT C K/ CPKEHF+RE ) )-PREVF) ANO 
Ado(DGAMF)<ABSCPREGAMF)) OR (HFQ > PREHFO ) THEN 
0 E G I N 

KLAGM* UHF«-( IF 'COUNTsl THEN 0 ELSE DHFO)* 
go TO L 5 END ELSE BEGIN F L A G «• 0 } 

FOK HO STEP 1 UNTIL ( a- A X I N C E X <- DL 0 M A X I NO E X ) DO 
BEGIN COMMENT IF THE ITERATION WAS UNSUCCESSFUL 

DISCARD IHE DATA AND ATTEMPT ONLE TO SATISFY 

terminal constraints; • 

PHI [ IUOL DPMI [ n ; 

LOUPL XPCL, 1 ]<-ULDXP[L, 13* 
end; 

coast+olocuast; 
ni^Ti+cOAsr; 
t f «■ r 2 + T B 2 ; 

IF I TER- 1 THEN begin 
writeupmge j ); 

COMMENT IF THE FIRST ITERATION IS UNSUCCESSFUL' THE 
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FOLLOWING MESSAGE 15 PRINTED OUT# 

W f-JI T E ( < " T H E TERMINAL CONSTRAINTS HAVE BEEN VIO"# 
"LAT ED EXCESSIVELY BECAUSE"#//# "THE NO"# 
"MlNAL CONTROL OR THE INITIAL CONDITIO"# 
"NS ARE"#//#"NOT CLOSE ENOUGH TO THE "# 
"OPTIMUM GUESS AGAIN">)#GO TO L ^ 

end; 

go to l 3 End; 

END 

ElSL BEGIN FLAG«-i; DHEMOOOOO/(2*EAIL)# COUNTS 
COMMENT THE FUlOrINu STATEMENT CONTROLS PROGRAM TERMINATION 
IF QHFSIOOO THEN PEGli\ 0 HP <- DHFO ; QU I T 4- 1 ; GO TO L5 

end; 



eno; 

L 5 * OVF«-SORT(</(XPLl#MAXlMTEX3<-DHF + RE))-XPC3#MAXlNDEXJI 

OGAMF«--XP£ A # MAX INDEX l; 

GO TO L 2 ; 

COMMENT Fl'»AL PRECISION RUN; 

L 6 1 A Jams ( A# 1# O# TF# 20. 20# P-A, 3 - 4 # XP# EUNCT) # 

COMMENT print 0U1 SENS I I I V I I I ES OF PAY-DEE TO ERRORS IN 
TERMINAL CONSTRAINTS# 

n r it e c i i i p n # 1 1 # r m p 1 1 # 2 1 ) ; 

COMMENT punch the OPTIMUM CONTROL HISTORY ON CARDS# 

HELL : NR1 TEC CARDS# H I 3TURY # FOR I <- q STEP 1 UNTIL MAX INDEX DO 

phi ud; 

GO 10 LIT la: END. 
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APPENDIX B. SECOND VARIATION COMPUTER PROGRAM 



begin comment e^nest c. lude»s box 219 second variation; 

REAL RE»GO/K,ISP,TF *T0VERWC»NUl*NU2*HH,V0*TEl*TE2* 

R* V,TSING*OLDPETM* DETM,0L0Te1'0LDTE2*HUK* T3,TA, 
OLDNIJ1* 0LDNU2, OLDT/ 

real SUM* FEE# M A SR AT/ S*C*Q,C1*C2*C3*C4* I NT j 
REAL ARR a v TEMP3#TEMP6ris31*TEMPA»TEMP5Ci:3*i:33# 

DELX, TEMPI* TEMP 2/ BC 1 : 3 ] * M I \ V , m* N C 1 : 3 * 1 : 3 3 * 
C0C0U*0!1153; 
integer ind* b k i n o ; 

INTEGER ITER*MAXINCEX# I*ONCE#K'<* j*l*st; 

INTEGER BOUND; 

REAL array OLDXP*XP*YPCOSA*0:i25]#PHI#HU*HUUtO»1253# 
rU/HUX c 1 j 3/ 0 : 1251 # ZZ[0:24,0:i25]* 
FL*SLCl*i#i;3*100:il03*ZPC0S3,0:n; 

LABEL L 1 * L 2 ; 

FORMAT HlST0RY(Afl8.12); 

FORMAT EICON THIS I T E R A T I ON"* // * X5* "N'J 1 = ", E20.ll* 
//* X5* 

" N U 2 = "*E20.li#//,X5*"TEl = "*E20,U*//* 

X 5 , " T E 2 = "* 

E20. 10,//,X5,"HUK = ", E20.ll); 

DEFINE LOOPl = FOR l«-l*2,3 DO *, LOOPJ = FOR J«-l,2,3 DO t» 
LOQPL = FOR L 4 - 1 , 2, 3 DO #; 



comment matrix inversion procedure goes here; 

PROCEDURE NOMINAL CT,X,OX); 

value t; real t; real array x*dxcu; 
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BEGIN 



Real s,c, masrat > tee# g; 

INTEGER I;Ji 

I «■ ENTIERCT/HH); J «■ I F I<MAxlNOEX THEN 1 + 1 ELSE Ii 
FEE <- CQ<- p Hl[I]) + (T/HH - I)x(PHl[J] - Q); 

R«-X[ 1 ]+re; 

masraTm-toverwoxT/isp; 

S*SINC XT A] )i C<-COS(XC A] ); V«-xC3i; 

DX [ 1 3 *• VxS; 

DX [ 2 3 «- VxC/R; 

D X [ 3 3 «■ ( q «- GoxTGVERwo/MA$RAT)xCOS(FEE) - KxS/R*2; 

D X [ A ] *- -K/V X C/P*2 + GxSIN(FEE)/V +DXC23; 

end nominal; 

PROCEDURE BACK ( T> Z'DZ )> 

value t; real t; real array z^dzcm; 

8EGIN 

REAL T9,kS>A,LAm1>LAM2>LAM3>KC>R2>R3/V2,S,C>SF>CF/ 

masrat, u> im; 

real array FX,hXX>F,Q, SS,M.*N>SM>FTN,FM>QN[ \ : 3 / 1 5 3 3 > C C , 
D/vtd,ntc/Hxuci ; 3i; 

INTEGER e; 

LABEL EN; 

TB<-TF-T; I«-ENT I EPCT6/HH+ .00000000 n; 

IF ISO THEN I <- c; MASRATM-TaVERWOxTB/ISp; 

J «■ IF kmaxindex then 1 + 1 else i; INT «-tb/hh - i; 
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R «- CU«-yp[ 1> I ] ) + INTxCXPt 1, J] - U) + Re; 
V «■ (U«-XPt3*I3) + IN1 x(XPC 3»J3 “ uj; 
s «■ SINC (U«-XPC4, 13 ) + INTx(xPU»J3 - U)); 
C «■ CQSC ( J«-XPC U, 13 ) + INTX(XPC«»J3 - U>>; 
SF <• SlN'( (U«-PHI t 13 ) + INTxcphICJD - U))i 
CF <- COS( (U“PHI [ I ] ) + INTxCPHICJ] - U))f 



COMMENT COMPUTE PARTIAL derivatives of system equations 
WITH RESPECT TO STATES/ 

FXU,13«-o; FX[1,23«-S; FX [ 1 , 3 ] <- VxC,* 

FXC2/1 ]*2xkxS/R* 3; FxC2,2D<-0; F X C 2 , 3 ] «— K xC /R *2 i 
FXC3/1]«--VxC/R*2 + 2xKxC/(VxR*3)J 

FX[3*23«-C/R + K*C/CVxR)*2-GOxTOVERhCxSF/CV*2xt/ASRAT> 
FXC3/3]«--VxS/R + KxS/(VxR*2); 

COMMENT COMPUTE PARTIAL DERIVATIVES oF SYSTEM EGUATlONS 

with respect tq control; 

FUC 1/ I 3*o; 

FUC2 / I ]«--( A<-GOxTOVERWO/MaSRaT)xSF; 

FUC 3/ I ]4-a/VxCF; 



comment compute derivatives of the hamiltonian; 



LAMUZC 223 ; L AM2*Z C 2 3 ] ; LAM3*Z[243; 
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HUC I J«-LAM2xFU[ 2* 1 3+IAM3XFUC 3* I3i 

HUUCI ] <--L a M2x AxCF-LAV3xAxSF/V - HUK; 

COMMENT CHECK HUU > O; 

IF HUUCI]>0 THEN 
IF T S I NG = 0 THEN 

begin 

TSI\G«-TB + HHx5; 

WRITEC<"HUU > 0 AT TIME = "/ F20 . 1 0>> TB ) ; 

write(<”huu = ’%e2c. io>/huuc n ); 
if tsing>tf then 

BEGIN 

tsing*o; 

TSING«-1/TSING 

end; 

GO TO EN 

ENO 

Else GO To en; 

HUXCl, I3«-HUXC3> I]«-HXUCl]*HXUC3]«-o; 
HXUC23«-HUX[2/I3«'-LAM3 xfU[3>I]/v; 
kS«-kxS; kokxc; p2«-rxr; v2«-vxv; r3«-r2xr; 
HXXCl>n«-»6xLAM2xKS/R2*2 + LAM3x(-6xKC/(R2xR2xV) + 2xVxC/R35; 
HXX[l,23«-HXX[?, 1 ] <--LAM3 x(2xKC/(V2xR3) + C/R2); 

HXX[ 1 # 3]«-HXXC3/l) <-2xLAM2xKC/R3+LAM3x(-2xKS/( R3xV > + VxS/R2 ) 
HXXC2,2]<--LAM3x(-2xAxSF/( V2xV)+2xKC/(R2xV2xV) )> 
HXXC2>33«-HXXC3»23«-LAMlxC + L AM 3 x ( -KS/ ( R2xV2 ) -S/R ) / 

HXXC 3# 3 ] *--l_AMlxVx$ + LAM2XKS/R2 + L AM 3x ( K C / C R 2x V ) - Vx C /R ) ; 



55 



kk«-o; 



LOOPL LOOP J 

begin 

FCL>J3«-FXCL*J3-FUCL*I3xhuX[J#I 3/HUUC I ]i 
QCL' J3«-FU t l, I] *FU[ J»I 3/HUUC n; 

SSCL*J3*HXX[L* J3-HXUCL3 xHXUCJ] /HUUCI3; 
M[L»J3«-Z[CKK«-KK + 1)]; N[l» J3«-Z[KK + 93i 
COMMENT STORE f AND SS MATRICES for USE In PROCEDURE NEWlAMDAJ 

if i >100 then 

BEGIN 

flcl, J> n-FCL> J]; ' 

SLtL# J/I3«-SSCL* J3i 
end; 
end; 

IF KMAXINDEX then 

DETM*-M[1,1]x(mC2,2]xmC3/3d-MC2,3]xmC3,23) + 
M[1,23X(M[2>33xM[3>13-Mc2,13xM[3>33)+ 
M[1*33x(mC2#13xhC3*23-m[2,23xmC3/13); 

comment check for change in sign op determinant of m; 

IF I<(maxINDEX-1 ) AND 

S1GN(0LDDETM)/SIGNCDETM) then 
IF T S I N G = 0 THEN BEGIN 

tsing^tb + H h y 5 ; 

WRITE(< m THE DETERMINANT of M changed sign at T = "# 
F20.10>,TB); GO TO EN 
END ELSE GO To EN; 
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A-HUt I3/HUUC I 3 J 
LOOP J BCG I N 

cccj]«-“Fucj# n x a ; 
o[J]<-hxucj] x a; 
end; 

olddetmoetm; 
loop j loo p l begin 
FMC J # L]*-0; 

ON C J,L3‘-0' 

SM[ J,L)*-0/ 
ftnc j,l3«-o; 

FOR KKM>2,3 DO BEGIN 

FM[J,L3«-FM[j,L]+ F[J,KK]x m [ k K / L ] / 
qn[j>l]«-qn'[j#l]+ qcJ/Kki* ncxk/Lj; 
rTN[J,L]«-FTN[J # L3+F[KK # J3xN[KK,L3/ 
SMCJ#L3«-SM[J,L3 + SSCJ#KK3 x M[KK,L3/ 

end; end; 

KK<-0; LOOPL MTDCL3«-NTCCL3«-OZCL*2l 3*-OJ 
LOOPL LOOP J BEGIN 

DZC (KK«-KK + 1 ) 3«--FM[L» J3+QNCL, Jj; 

DZCkk + 93 4 " SMCl>J3 + FTNCL'J}; 

DZC CC«-L + 2l ) 3<-DZCE3 + FX C J, L 3 * Z I J + 2 1 3 ; 

MTD[L3 < -MTDCL3 + M[J,L3XDCJ3; 

NTCCL3«-NTC[L3 + NCJ>L3 X CCCJ3> 
end; 

loopl CZCL + lB3«--'MTDfLl + N T C C L 3 ; 

en: end back; 

procedure controlc t» x>dx ) ; 



57 



value t; real t; 

BEGIN 



real array x>dxcij; 



Lis 

L2S 



COMMENT 



LABEL L1,L2>L3,LR; 

IE T>ONCE x . 9999999 THEN 
BEGIN 

IE BOU\'D = 0 THEN 

begin 

IND*-ENTIER(ONCE/HH + .0000001 ) J 
GO TO L2 

end; 

i nd«-entier(once/hh+. oooooo n + i; 
once«-once+hh; 

B<lNn«-MAXlNDEX-lNDJ IF InD>MAXINDEX - 7 THEN GO TO LA 

k k <- o ; 

loopi begin 

B[ il^ZZCI+lS/BKIND]; 

LOO 3 J BEGIN 

MINV[I»J]«-MCl>J]«-ZZ[(KK«-KK+i)#BKlN0 3j 
N[ 1/ J3«-ZZCKK + 9»BKIND3 

end; 

end; 

invertc minv> 3/ 1 ); 

IE 1=1 THEN WrITE(<”M IS SINGULAR AT T = ",E6.2>/T); 
COMPUTE INITIAL CONDITIONS FOR PROCEDURE NEWLAMDA; 

IE I ND = 1 00 THEN 
LOOPI BEGIN 

L 0 0 P J BEGIN 
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COMMENT 



TEMPfltl, J3«-MINV[I, J3J 
TEMP5CI/J3*- N[ I, J3; 

END/ 

TEMP6C I 3 «- B C I 3; 

end; 

IF I N D = 1 0 2 THEN 
begin 

DELXC13«-XP[1>1003“OLOXP[1»1003; 

DELX[23«-XP[ 3 / 100 3-0UDXPC 3* 100 3; 
DElX[33«-XP[4/1003-0LDXPU/1003/ 

LOOP I BEGIN 
5UM«-o; 

LOQPJ 5UV«.SUMrTEMp5C J# 1 3xDELX C J3 i 
TEMPI [ I 3«>SUM*TEMP6t I 3; 

END/ 

LOOPI BEGIN 

sum*o; 

LOOP J SUM«-SL'M + TEMpac J, I 3><TEMP1[ J]; 

zpc i/ 03 «-sum; 

END/ 

end; 

compute COEFFICIENTS FOR COMPUTING new control; 
LOOPI BEGIN 

sum*o; 

LOOPJ SUM^SUM+FUC J/ IND]XMINVC 1/ J3/ 

TEMP 1 [ I 3 «• SUM / 

end; 

LOOPI BEGIN 



59 



sum«-o; 

LOOP j SUM4-SUK + TEMP1C JJXK'C j]; 

TEMP2C I 3*SUM; 

end; 

L«* IF INO >MAXINOEX THEN GO TO L3; 

TEKP3C 1 3 «-OLDXp[ 1 , INOl J 
TEMP3t2]«-OLDXPC 3, INDi; 

TEMP3C 33<-3LDXPt A, IN03; 

T3«-T4«-o; 

LOOP J BEGIN 

COC J, I\'D3«--(HUXt J / IN03 + TEHP2C J3 3/HUUt I NO 35 
T3«-T3 + TEMP1 1 J3XZZC J + 18 /BKINd3; 

end; 

T3«-(T3 + HUtIND3)/HUUCIND3i 

LOOP J T4«-T4+C0[J> IND]x3EMP3[J]; 

C0tA»IN03«'PHItIND3-TA-T3; 

WRITEC<"INDEX = "»IA>#IND3J 

WRITE(<"C1 = H /E20. 10»"C2 = " , E 20 . 1 0 . "C 3 = ",E20.10> 
= "/E 20 . 10 >, FOR 1*1 STEP 1 UNTIL A DO cOCI»IND33 
WRITE! T3/ FU 3/ 

IF B 0 U N D = 0 THEN 

begin 

B0Un~i«-80UND + 1; 

GO TO LI 
end; 
end; 

L 3 : I «- ENTIERCT/HH3; J «■ I F I<MAxlNDEX THEN 1 + 1 ELSE I 

IF T/OLDT THEN 
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BEGIN 



oldt«-t; 

int«-t/hh-i; 

Cl«-(Q«-C 0 Cl/I 3 ) + lNTx(C0tl/J3-Q); 

C2«-(0«-C0C2/I]) + INTx(c0C2/J]-0); 

C3«-(Q«-COC3,I]) + lNTx(COt3>J3-Q)i 
C4 *-cq«-cou,i3) + initx(cou,j]-g); 
end; 

FEE«-pHl [T/HH3«-ClxXr 13+C2xX[33+c3xX[43+C4; 

R«-X[ 1 ] + re; MASRATM-TOVERWOxT/ISP; 

S«-SIN(X[43)i C«-C0SCX[/4]); V «■ X C 3 3 ; 

oxen*- vxs; 

DXC23*- VxC/R; 

D X C 3 3*- ( 3 <-G 0 xTCvERW 0 /MASRAT)xCQSCFEE) - Kx£/R*2/ 
DXU3«--K/Vxc/R*2 + GxSIN(FEE)/V + D X C 2 3 J 

eno control; 

COGENT integration procedure goes here; 

PROCEDURE NEWLAMDA(T,LAM/OLAM); 

value t; real t; real array lam/Dlamcu; 
begin 

real a»int,q; 

real ARRAY D/DELX/SQX >FTDL[ 1 * 33,SDELX[ 1 : 3, 100: 1 10] 
/FlICI *3/1 *33; 

integer ind/.kk/Jj; 

LABEL L1/L2*L3; 

IF T >0 NC E x , 999 Q 90999 THEN 

begin 
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IF B0UND=0 THEN 
BEGIN 

i nd«-entier(once/hh + . ooooooi >; 

GO TO L 2 

end; 

L 1 1 IND*-ENTIER(0NCE/HH+.0000001 ) + i; 

« 

once«-once+hh; if ind>haxindex then go to l3; 

L 2 S DELXC1]*-XP[WINd3-0L0XPC1>IND3; 

DELX[23«-XP[3, IN0 3-0L0XPC 3» I NO 3/ 

OELXC 3 3<-XP[ H, INO 3-OLD XP C I NO]; 

A*HUC IN03/HUUCIND3; 

L 0 0 p I BEGIN 

or 1 3 «-hux c i> ind3xa; 

SOELXC I, IMD3«-o; 

LOOP J SOELXC l> lND3«-SDELXt I > IND3+SLC l> J> IN03X 
delxc ji; 

SDELXC I , IfJO j«-50ELX[ I > IND3-D[ 1 3> 

end; 

IF 30Un0=0 then 
BEGIN 

BOUND«-BOUNO + i; 

GO TO U 

end; 

end; 

l 3 1 kk«-ent ierc t/hk ) ; jj<-if kk<maxindex then kk + i else kk; 

int«-t/hh-kk; 

LOOP I BEGIN 

SDX[ I3«-(Q«-SDELX[I>KK] ) + INTxCSDELXr I* JJ3-Q); 
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loop J flici,J3«-(q«-fi[I/J,kk]) + intx(flci*j>jj3-q)J 
end; 

LOOP I BEGIN 
FTDLC I 3«-o; 

LOO p J FTDL[I3«-FTDL[I3+FLI[J>I3 *LAM[J]; 

L'LAMt I 3«- -SDXC n-FTDLC I ]; 
end; 

end newlamda; 

COMMENT MAIN PROGRAM BEGINS HERE; 

RE«-2.09G7; GO«-3?.17; k<-g0*re*2; iTer*o; I Sp*-300; 

0LDT«-6900; 

READ(XFC3 / 03#XPt6,03»TF»T0VERW0*Nui*NU2/HH,HUK); 

XPCl>C3«-XP[2*03«-0; 

XP[4,03«-Xp[ft, 03/ 57. 2957795131; 
m;.XINDEX<*EMTIERCTF/HH + ,51); 

READC hi STDRY» for I +0 STEP 1 UNTIL maxinDEX DO PHICI3); 
WR1TECHISTqry> FOR I*0 STEP 1 UNTIL MAXINDEX DO PHICI3); 

comment integrate system equations with nominal control; 

ADAMSU/HH, 0> T f » HH» TF/ 0 , 0/ XP, NOMINAL); 

vo<-sqrt(K/(x p ci>maxindex]+re)); 

N02«-NU2/C2xV0>; 

TE1«--XPU, MAXINDEX3; 

TE2«--V0 +XPC 3/MAXINDEX3 i 

L 1 : IF I TER/ 0 AND ABSCOlDTEI )<ASS(TE1 ) AND ABS ( OLDT E2 ) < A3S ( TE2 ) 
AND OLDXPC l.v,AXlNOEy]>XPC 1 * MAXINDEX3 THEN 
BEGIN 

TEUTEl/2; TE2<-TE2/2; HUK<-HUKxlQ; 
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END ELSE HUK«-HUK/2; 

COMMENT COMPUTE INITIAL CONDITIONS FOR BACKWARD INTEGRATION 
r«-xpci,maxindex3 + re; v«-xpc 3> max index]; 
for i«-2 stfp i until ia do zzci/03*o; 

ZZC l/03*-2xVi 
Z ZC4/0] <- k/R*2; 

ZZCIO^O] •- -AxVxKxNU2/R* 3 ; 

ZZC12'03 «- -ZZC4>0i; 

ZZC13/0] « -2xKxNU2/R*2; 

ZZC 15>0j <- + ZZ[l/03; 

ZZC17>0] «- 15 
ZZC19/0]* -TEi; 

ZZ [ 20> 0 3 «- -TE2; 

ZZC21/03*0; 

IF I TER=0 THEN BEGIN 
ZZt22»03*j+K>«NU2/R*2i 
ZZ L 2 3' 0 3 «■ 2 xVxNU 25 
ZZ t 24 > 0 3 * "NU 1 / 

END ELSE BEGIN 

bouno*o; 

DNCE*200; 

A DAM SC 3,HH, 2 00' TF, TF, JF , P-4, ZP > N'EWL AMD A ) ; 

ZZ[22*03*ZPtl> 13+ZZC22/ 03; 

ZZC23'0]«-ZP[2' 13+ZZC23»03J 
ZZt2A,.03<-7.P[3>n + ZZ[24,035 

end; 

TSI NG<-0; ONCE*0; 

detm*o; 
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convent integrate m,n,b, and lambda equations backward; 

ADAMSC2A* HH> 0, T F > HH» 10, 0 /O > ZZ>BACK); 

CONVENT STORE OLD VALUES QF STATES BEFORE RUNNING NEW TRAJECTORY; 
FOR I<-0 STEP 1 UNTIL WAXINDEX DO 
FOR L<-1>2,.3,4 DO OLQXP C L > I 3 «-XP t L M 3 » 

COMMENT IF DET(m) changed sign then A CONJUGATE point EXISTS 
or IF Huu WENT POSITIVE DURING THE INTEGRATION THE LEGEnDRE 
CONDITION IS NOT SATISFIED, TSINg IS SET TO The TIME at 

which either occurred and the porward integration is started 

AT THIS POINT rather THAN ZERO-* 

0NCE<-0; 

bound«-o; 

IF T S I N G = 0 THEN AD Ay, S t a » HH , o , TP » HH * TF , 0 > 0 * XP , CON TR OL ) 

ELSE BEGIN 

if tsing>tf then begin 

WRITEC<"HUU WENT POSITIVE OR A CONJUGAL 
"TE POINT OCCURRED TOO CLOSE TO THE > 

"END or THE TRAJEcTGHY">);GO to L2 end; 
ONCE^-ENT IERC TSING+. 1); L<-ENTIER( TSING/HH+, 1 ) J 
FOR J«-l>2»3>4 DO YP[J>0]*-XP[J,Lj; 

ADAM S(4.HH,TS!NG/TF,HH>HH,0,C,YP, CONTROL); 

FOR I-L STEP 1 UNTIL MAXINQEX DO 
FOR J<-1,2,3,4 DO XP[J#I]«-YP[J/I-L3i 
END^ 

WRITECFOR I<-0 STEP 1 UnTIl MAXInDEX DO PHICIJ); 

VO«-SQRT(K/(XPt UMAX INDEX ] + RE) )J 
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OLDTElMEl ; OLDTE2«-TF.2; 

TEl«--XPLfl,MAXlNOrx]; 

TE2<--V0 +XPC3*MAXlNDEX] ; 

comment if improvement in terminal errors and p ayoff is 

LESS than EPSILON ( USERS CHOICE) THEN STOP THE PROGRAM; 

IF A35C0L0TE1-TE1 )<.OCl and A 3 S ( OLD T E 2 - T E 2 ) < 1 

AND ABS(OLDX p n/MAXJNCEX]-XPll»MAxlNDEX] )<500 THEN 

BEGIN 

WRITECFOR 1*1 STpP 1 UNTIL MaXiNDEX DO PHIII3)) 

GO TO L 2 
END ELSE 
BEGIN 

ITERMTER + 1* 

WRITE(F1>NU1'NU2>TE1,# TE2>HU 1 <); 

GO TO LI 

end; 

L2: end. 
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